/*matrix.cpp*/
#include "matrix.h"

//设置大小,分配空间
Matrix::Matrix(int row, int col) {
	Row = row;
	Col = col;
	matrix.resize(row);
	for (int i = 0; i < row; i++) {
		matrix[i].resize(col);
	}
}
Matrix::Matrix(vector<double> v,bool div) {
	int n = v.size();
	if (div) { this->Col = 1; this->Row = n; }
	else { this->Col = n; this->Row = 1; }
	this->matrix.resize(this->Row);
	for (int i = 0; i < this->Row; i++) {
		matrix[i].resize(this->Col);
		for (int j = 0; j < this->Col; j++) {
			matrix[i][j] = v[i * this->Col + j];
		}
	}
}
void Matrix::resize(int row, int col) {
	this->Row = row;
	this->Col = col;
	this->matrix.resize(row);
	for (int i = 0; i < row; i++) this->matrix[i].resize(col);
}
//获得矩阵行数
int Matrix::row() {
	return this->Row;
}
//获得矩阵列数
int Matrix::col() {
	return this->Col;
}
//区块设置数据
void Matrix::setblock(int rsta, int csta, int rsize, int csize, double* pBlock) {
	for (int i = 0; i < rsize; i++) {
		for (int j = 0; j < csize; j++) {
			matrix[rsta + i - 1][csta + j - 1] = pBlock[i * csize + j];
		}
	}
}
void Matrix::printm(int n) {
	for (int i = 0; i < this->Row; i++) {
		for (int j = 0; j < this->Col; j++) {
			std::cout << std::setw(n) << this->matrix[i][j] << " ";
		}
		std::cout << std::endl;
	}
}
//矩阵乘法
Matrix Matrix::operator*(const Matrix& mrig) {
	if (this->Col != mrig.Row) exit(1);
	Matrix mequ(this->Row, mrig.Col);
	for (int i = 0; i < this->Row; i++) {
		for (int j = 0; j < mrig.Col; j++) {
			for (int k = 0; k < this->Col; k++) {
				mequ.matrix[i][j] += this->matrix[i][k] * mrig.matrix[k][j];
			}
		}
	}
	return mequ;
}
Matrix Matrix::operator*(const double& num) {
	Matrix mequ(this->Row, this->Col);
	for (int i = 0; i < this->Row; i++) {
		for (int j = 0; j < this->Col; j++) {
			mequ[i][j] = this->matrix[i][j] * num;
		}
	}
	return mequ;
}
//正定对称矩阵LDL^T分解,储存为[L U]
Matrix Matrix::LDLdecom() {
	if (this->Row != this->Col) exit(1);
	int N = this->Row;
	Matrix mdecom(N, N);
	double sum = 0;
	for (int n = 0; n < N; n++) {
		mdecom.matrix[n][n] = this->matrix[n][n];
		if (n == 0) sum = 0;
		else {
			for (int i = 0; i < n; i++) sum += mdecom.matrix[n][i] * mdecom.matrix[i][n];
		}
		mdecom.matrix[n][n] -= sum;
		sum = 0;
		for (int k = n + 1; k < N; k++) {
			mdecom.matrix[n][k] = this->matrix[n][k];
			if (n == 0) sum = 0;
			else {
				for (int i = 0; i < n; i++) sum += mdecom.matrix[n][i] * mdecom.matrix[i][k];
			}
			mdecom.matrix[n][k] -= sum;
			sum = 0;
			mdecom.matrix[k][n] = mdecom.matrix[n][k] / mdecom.matrix[n][n];
		}
	}
	return mdecom;
}
//系数为正定对称矩阵的方程求根
void Matrix::root(vector<double>pb, vector<double>&px) {
	int N = this->Col;
	if (N != this->Row) exit(1);
	Matrix mldl = this->LDLdecom();
	double sum = 0;
	for (int i = 0; i < N; i++) {
		px[i] = pb[i];
		if (i == 0) sum = 0;
		else {
			for (int j = 0; j < i; j++) sum += mldl.matrix[i][j] * px[j];
		}
		px[i] -= sum;
		sum = 0;
	}
	for (int i = N - 1; i > -1; i--) {
		if (i == N - 1) sum = 0;
		else {
			for (int j = i + 1; j < N; j++) sum += mldl.matrix[i][j] * px[j];
		}
		px[i] = (px[i] - sum) / mldl.matrix[i][i];
		sum = 0;
	}
}